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'"^p 1 A central question in developmental biology is how cells interact to organize into 

tissues? In this paper, we study the role of mesenchyme-ectoderm interaction in the 
growing chick limb bud using Glazier and Graner's Cellular Potts Model, a grid- 
based stochastic framework designed to simulate cell interactions and movement. 
We simulate cellular mechanisms including cell adhesion, growth, and division and 
diffusion of morphogens, to show that differential adhesion between the cells, diffu- 
sion of growth factors through the extracellular matrix, and the elastic properties 
of the apical ectodermal ridge together can produce the proper shape of the limb 
\Q \ bud. 
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1 Introduction and biology 



The chick Gallus gallus is a major model organism in developmental biology. 
As a common and accessible source of embryos, chicks attracted the interest 
of the ancient Egyptians, as well as of Hippocrates and Aristotle. Aristotle's 
studies influenced embryo logists until the 17th century, including Albertus 
Magnus, Harvey, Malpighi, von Haller, and Wolff [1]. Embryology progressed 
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in the 19th century, partly thanks to technical advances applied to the chick 
embryo (Pander and von Baer, followed by von Ebner, Hensen, Rauber, Koller, 
and Remak). Modern developmental biology was born at the beginning of the 
20th century, when experimental embryology combined with developmental 
genetics (Roux and Morgan). The chick remained important to experimental 
embryology at this time (Rawles, Fell, Rudnick, Graper, Wetzel, Adelmann, 
Pannett, Huxley, De Beer, Weiss, and Waddington). From the 1940s to the 
1970s, studies of chick development lost momentum as attention shifted to 
model organisms such as Xenopus laevis which were better suited to early 
molecular biology methods. New technologies, including in vivo electropora- 
tion, embryonic stem cell separation, transgenesis, and the completion of the 
chick genome sequence, have restored the appeal of this versatile experimental 
organism [2]. 

Limb development is an excellent model of vertebrate embryo patterning [3]. 
Chick embryos permit manipulation of the developing limb in vitro or through 
a window in the egg shell in vivo. Saunders' two key experiments on limb 
development in the chick showed that the removal of the apical ectodermal 
ridge (AER), a narrow band of specialized ectodermal cells located at the tip 
of the growing bud [4], caused distal truncations in the limb [5], and that 
the grafting of the posterior zone of polarizing activity (ZPA) to the anterior 
margin led to the development of additional digits [6]. These experiments 
defined the cell interactions in the developing limb that control growth and 
patterning along each of the three limb axes, and led to the formulation of 
various models of these interactions [7,8,9]. 

Hamburger and Hamilton published a detailed atlas of the stages of chick de- 
velopment [10]. The mesenchymal tissue of the early limb bud forms a paddle- 
shaped mesoblast with a surrounding layer of epithelium (the ectoderm), as 
shown in Fig. 1. The interaction of the distal portion of the mesoblast (the 
apical zone) with the AER drives proximodistal development of the chick 
limb bud. The AER secretes diffusing fibroblast growth factors (FGF) [11], 
which keep the apical- zone mesenchyme in a noncondensed and proliferating 
state [12,13], control cell division [14,15], and shape the bud [9]. 

Formation of the limb bud starts when the lateral plate mesoderm cells in 
the animal's flank begin to proliferate, approximately 1.5 days after fertiliza- 
tion [16,17,18]. Mesenchymal cells secrete FGF10, the key mesenchymal factor 
in limb bud initiation, which induces formation of the AER [19,20]. The AER 
then sustains the outgrowth of the limb by secreting FGF8 [21]. Applica- 
tion of FGF8 to the flank of young chick embryos induces the development 
of additional limbs [21]. Moreover, FGF8 can replace the AER to maintain 
limb growth, while misexpression of FGF8 causes severe limb truncations [21]. 
FGF8 also maintains secretion of FGF10 in the mesenchyme, and positive 
feedback between FGF8 and FGF10 is essential for continuous growth of the 
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Fig. 1. (a) The chick embryo at Hamburger and Hamilton stage 18. The rectangle 
encloses the area shown in panel (b). The embryo is approximately 7-8mm long, 
and the width of the limb bud (~ 1mm) remains almost constant as the limb 
grows. Adapted, with changes, from Ref. [29]. (b) A larger picture of the chick wing 
bud at Hamburger and Hamilton stage 16, showing somite numbers, mesoderm, and 
ectoderm. Adapted, with changes, from Ref. [29]. (c) A sketch of the embryonic limb 
bud showing the major axes, the apical zone (light grey), precartilage condensations 
(medium grey), and cartilage (dark grey). Adapted, with changes, from Ref. [30]. 

chick limb [22]. Thus, FGF8 appears to be the key molecule in outgrowth 
and patterning of the developing vertebrate limb [21]. Applying beads soaked 
in FGF1, FGF2 and FGF4 to the flank also induces additional limbs [23]. 
However, biochemical and expression studies suggest that none of these is the 
endogenous signal that initiates limb development [21]. 

4 days after fertilization, the mesenchymal cells condense and begin to form 
cartilage. Skeletogenic condensation is a multistep process that requires ini- 
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Fig. 2. (a) Expression of FGF8 (black) in the AER of a Hamburger and Hamilton 
stage 21 chick wing bud. Adapted, with changes, from Ref. [21]. (b) Schematic 
representation of the interaction between the AER, ZPA, and mesenchyme in the 
early chick limb. 

tiation, establishment of boundary conditions, cell adhesion, proliferation, 
growth, and cessation of growth (for a review, see Ref. [24]). Extracellular 
matrix (ECM) macromolecules, such as fibronectin, tenascin and transform- 
ing growth factor (3 (TGF/3), accumulate at sites of precartilage condensation 
changing cell-ECM interaction [25,26]. To explain skeletogenesis we must find 
a dynamical model of the spatio-temporal distribution of these molecules that 
can generate the skeleton [27]. If skeletogenesis involves a Turing instability of 
reacting and diffusing chemicals [28], we must model its boundary conditions 
(the dynamics of limb shaping, cell division, and cell growth), which control 
where and how the cells aggregate [27,30,31]. 

The ZPA regulates anteroposterior limb patterning by secreting sonic hedge- 
hog (Shh) [32]. Shh initiates and activates expression of secondary signaling 
molecules, including bone morphogenetic protein 2 (BMP2) in the mesoderm 
and FGF4 in the posterior part of the AER [33]. BMP2 mediates long-range 
effects of Shh, while FGF4 maintains the production of Shh in the ZPA. Fig. 2 
shows the expression of FGF8 in the AER, and the schematic interaction be- 
tween the mesenchyme, AER, and ZPA. We know less about patterning across 
the dorsoventral axis. Experiments suggest a role for ectodermal signalling in- 
volving Wnt7a [34] which also helps maintain Shh expression in posterior 
mesenchyme [35]. 

In this paper we concentrate on the shape of the growing chick limb bud. We 
assume that the AER secretes FGF8 (which decays with a given life-time in 
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mesenchymal tissue) at a constant rate, to account for the regulatory role of 
FGF10. FGF8 diffuses into the underlying mesenchyme and induces the cells 
to grow and divide. We omit the role of the ZPA and condensation, leaving the 
anterior-posterior asymmetry, as well as cartilage formation, to future work. 
For simplicity, we also assume that the limb bud is dorsoventral symmetric. We 
show that cell adhesion and growth together with diffusion and decay of FGF8 
can generate the proper shape of the bud, if we choose the physical parameters 
properly. We also show that the AER is crucial not only for maintaining the 
outgrowth of the limb (a permissive function), but also for forming the desired 
shape (an instructive function). 



2 Computational model 

Our simulations use CompuCell3D [36,37,38,39,40,41], a multi-model software 
framework for simulation of the development of multicellular organisms ( mor- 
phogenesis ). CompuCell3D models the interaction of the gene regulatory net- 
work with cellular mechanisms such as cell adhesion, growth, and division, 
diffusion of morphogens, chemotaxis (movement of cells along a gradient of 
a chemical diffusing in the ECM), and haptotaxis (movement of cells due to 
changes in cell-cell and cell-ECM adhesion which depend on the concentration 
of cell adhesion molecules). CompuCell3D implements Glazier and Graner's 
Cellular Potts Model ( CPM) [42] , a grid-based, stochastic model, which de- 
scribes biological mechanisms in terms of effective interactions and constraint 
energies [43] . Some recent applications of the CPM include studies of chicken 
retinal cells [44], Dictyostelium discoideum [45], cancer invasion [46], flowing 
foams [47,48], and viscous flow [49]. 

iV spatially distributed cells, indexed by a, lie on a lattice, and the value at a 
lattice site (pixel) (i, j, k) is a if this site lies in cell a. A collection of connected 
lattice sites with the same index represents a cell. The evolution of the cells fol- 
lows modified MonteCarlo-Boltzmann-Metropolis dynamics [42,43], with the 
cells rearranging their positions to minimize their total free energy [50]. 

In our simulations, the CPM Hamiltonian includes three main terms: 

(i,j,k),(i>,j>,k>) 

+ E MrW(a) - V t (a) f + £ \ s (r)(S(a) - S t (a)) 2 , (1) 

a a 

where r denotes the type of a cell. The first term describes the surface adhesion 
energy between cells and between cells and their environment, and the surface 
tension coefficients are symmetric, J T y = J r >, T - The second term corresponds 
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to the compressibility of the cells: V is the actual volume of a cell and V t its 
target volume. We simulate growth of the cells by increasing their V t . The 
third term represents the elasticity of the cell membrane: S is the surface area 
of a cell and S t its target surface area. In a two-dimensional CPM we omit V 
and V t , as in Ref. [42,43]. The volume constraint prevents cell disappearance, 
and the surface term fixes the average shape of the cells. We use the simplest 
form of surface term [36,37], although a term of form (SV~ 2 / 3 — StV t 2 ^ 3 ) 2 is 
more consistent with classical elasticity [51]. 

At each step we randomly select a lattice site (i, j, k) and change its value from 
a to the index of an arbitrary lattice site from the first-order neighborhood (a 
lattice site which is connected to a given pixel by a side) a', with a probability 
P: 



where AH is the difference in energy produced by the change, 6 is the Heavi- 
side step function, and T is a parameter ("temperature") corresponding to the 
amplitude of cell membrane fluctuations [52]. To calculate the energy result- 
ing from cell-cell adhesion interactions we consider first and second nearest 
neighbors, which reduces lattice anisotropy effects compared to first nearest- 
neighbor calculations. One Monte Carlo Step (MCS) corresponds to n flip 
attempts, where n is the total number of lattice sites. 

Our simplified limb bud includes 5 types of cells: mesenchyme (M), AER (A), 
regular ectoderm (e), regular mesoderm (m), and the medium (w) in which 
the bud is immersed. The AER secrets a growth factor (FGF8) which diffuses 
into the mesenchyme and mesoderm. The secretion occurs at a constant rate 
s, in any pixel containing the center of mass of an AER cell. The constancy of 
s ensures that the length of the apical zone, where the mesenchymal cells grow 
and divide, remains approximately constant [53]. We implement diffusion on 
the CPM lattice using the forward Euler method at each MCS. We use no-flux 
boundary conditions at the walls of the box which contains the growing limb, 
and at the boundary between the limb and the medium. The mesenchymal 
cells respond to the FGF8 concentration and grow at a rate proportional to 
this concentration, 



where a is a constant. We do not impose a priori any growth on the AER, 
ectodermal, and mesodermal cells, 



P(*(i,j, k)^a{ 



i'j f , k')) = 8{AH)e- AH/kT + e(-AH) 



(2) 



lV) = a[FGF8], 



(3) 



V t (A) = Vt(e) — Vt(m) = 0. 



(4) 
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In fact, the AER cells do not grow [9,54]. Our simulations show (see the next 
section) that large Jm w induces the ectodermal cells to grow and divide so 
that a constant monolayer remains on the surface of the bud. 

Fig. 3(a) shows the initial (3D) configuration. The dimensions of the box 
(in pixels) are 275x200x25 (the z— axis being perpendicular to the page). 
For simplicity, this configuration is dorsoventral symmetric, which means that 
the AER is represented by a thin plate rather than a stripe. Each cell initially 
occupies a cube 5x5x5 pixels, and 1 pixel approximately corresponds to 2/rni. 
A cell which doubles its volume with respect to its initial value, i.e. reaches 
V = 250, undergoes mitosis and splits into two cells with equal volumes. 
We estimate the diffusion constant for FGF8 moving in the ECM from the 
Einstein formula, 



where k B is the Boltzmann constant, T a is the absolute temperature, and kd rag 
is the drag coefficient in the relation between the drag force acting on a body 
moving slowly in a viscous fluid and the velocity of the body, F visc = k drag v. 
For a spherical body kd rag = 6nr]R, where rj is the viscosity of the fluid, and 
R is the radius of the body. 

For FGF8, which is an average-size protein, we have R ~ 10 _9 m. In limb 
bud ECM we have rj ~ lO 2 ^ [55,56]. Therefore, at room temperature we 
estimate D FGFS ~ 2-10~ 15! y-. The diffusion constant sets the time scale of our 
simulations, which measure time in Monte Carlo Steps, and distances in pixels. 
We set 1 MCS ~1 min, which yields D ~ O.l 2 ^^. Experiment shows that 
the proximodistal size of the apical zone L is approximately 200/im [53]. Our 
simulations set this length by appropriately chosing the FGF8 decay constant, 
kFGFs- From the diffusion equation we obtain the following relation: 



which also results from a simple dimensional analysis [30]. In the simulation 
shown in Fig. 3, we chose k = 0.001 MCS -1 , which corresponds to L ~ 50^m. 
Rescaling the length of the apical zone to the value 4 times smaller than 
the experimental one greatly speeds up our simulations. However, it should 
also produce the limb bud 4 times smaller (the cells keep their real size). 
This rescaling changes the relation between the CPM parameters and the 
corresponding physical values, if their units contain time (such as viscosity), 
since the time scale is set by the diffusion constant, but does not change the 
dynamics of limb growth. 



D = 



k B T a 



drag 




(6) 
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3 Results and discussion 



All parameters appearing in formulae (l)-(3) have CPM units (pixels, MCS). 
The "temperature" T sets the values for the surface adhesion coefficients J TT i 
and the volume constraint coefficient Ay. Here, we choose these values to 
prevent cells from disappearing or freezing, and use the hierarchy of surface 
adhesion coefficients to keep the cells from desorting. Thus, for each pair of cell 
types r and r' we set J TT < J TT / and J T / r / < J TT /. We remind the reader that, 
in the CPM, stronger adhesion corresponds to a lower value of J. Moreover, 
Jaa must be very low so that the AER does not break. The same is true for 
J ee , although J ee can be larger than Jaa since the ectoderm stretches during 
limb outgrowth. Since embryonic tissues are viscoelastic materials, we could, 
in principle, obtain their surface adhesion coefficients from parallel plate com- 
pression experiment [57,58,59,60]. Moreover, we could find the shear modulus 
which sets the value of As, and viscosity. Of the surface tension coefficients, 
we know so far only J mm [56,59]. The values of s and a are not independent, 
since their products drives the growth of the mesodermal cells. Thus, without 
loss of generality we may set s = 1, and we choose a = 4 so that the interval 
between cell divisions is of the order of lOh which is consistent with the actual 
duration of the chick limb mesenchymal cell cycle, which is 13 — 14h [61]. 

We use the following values for the cell-cell and cell-medium adhesion coeffi- 
cients, which satisfy the above inequalities: 



J ww — 0, Jaa — 1, Jmm — 4, Jam — 20, Jaw — 12, Jm w — 120, 

Jew = 12, J Ae = 12, J Me = 20, J ee = 6, J m w = 40, J Am = 30, 

Jmtti — 10, Jem = 20, J mm = 6. (7) 

They are of the same order of magnitude as those in Ref. [62]. We also set 
T = 30 following the same reference. Our choice of parameters is somewhat 
arbitrary and only experimentation can provide their exact relative values. 
The initial value of V t for all cells is equal to the initial volume, 125, and we 
keep S t V t ~ 2 ^ 3 = 5 constant, so that the cells remain round and do not stiffen 
as they grow. We chose Xs = 0.1 which suffices to round the cells. For all types 
Ay = 5, except for the AER which has Ay = 10. For these parameters, the 
AER is less elastic than other tissues, preventing the ridge from stretching and 
increasing the size of the apical zone, which would result in a bulbous bud [9]. 
We treat the medium as a large cell without volume or surface constraints. 

Fig. 3 shows the results of the corresponding simulation. The width initially 
grows rapidly, but the growth decreases and the width remains approximately 
constant. At Hamburger and Hamilton stage 17, which the chick passes be- 
tween 52h and 64h after fertilization, both wing and leg buds lift off the 
blastoderm. This moment corresponds to MCS in our simulation. At stage 
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23, after approximately 4 days, both buds are approximately as long as they 
are wide. In our simulation, this stage corresponds to ~5000 MCS and the 
width is approximately 100 pixels which is in agreement with the expected 
value of 200/xm. Thus, 5000 MCS correspond to 1.5 days, giving the relation 1 
MCS ~0.5 min which is consistent with the relation we obtained by setting the 
CPM value of the diffusion constant (1 MCS ~1 min). In Fig. 4 we show the 
configuration at 1000 MCS (from a different simulation) and the correspond- 
ing FGF8 concentration distribution, to demonstrate the underlying factors 
shaping the bud. 

The most important parameters in our simulation are: the rate of FGF8 se- 
cretion, its diffusion constant, and its decay constant. The first two cannot be 
too small, nor the last one too large, otherwise the FGF8 concentration will 
not suffice to make the mesenchyme cells grow. However, if one of the first two 
parameters is too large, or the decay constant is too small, the anteroposte- 
rior size of the apical zone increases, resulting in a bulbous bud, as in Fig. 5. 
The AER cells must not grow or divide, otherwise, even for slow growth, the 
size of the apical zone would grow due to the increasing number of AER cells 
secreting FGF8 (Fig. 6). For the same reason, the values of Ay and X s for the 
AER cannot be too small. We repeated this simulation for different values of 
Jik, Ay, and As, obtaining similar shapes of bud. 

Fig. 3 also shows the effects of lattice anisotropy on the shape of the limb. 
The tip of the bud initially aligns with the AER at an angle of 45°, forming a 
trapezoid. In time, the angle of alignment increases since only the most distal 
cells grow. Eventually, the tip becomes a polygon with more sides, while the 
rest of the bud remains smooth. This anisotropy appears when we limit spin 
flips in the CPM to nearest neighbors. To resolve this problem, we could extend 
spin flips to the second-order and third-order neighborhood (lattice sites which 
are connected to a given pixel by an edge and a corner, respectively). 

In Fig. 7 we present the time-dependence of the length and width of the limb, 
and the length of the AER. The length of the bud is approximately a lin- 
ear function of time, which is consistent with our assumption of a constant 
secretion of FGF8 and of mesenchymal growth proportional to the FGF8 con- 
centration. We measure the width of the bud as the maximal width of the 
mesenchyme. It increases most rapidly at the beginning. Later, growth slows 
to zero since the proximal mesenchyme ceases to grow. The length of the AER 
behaves similarly, except it grows fastest growth a little later. 

Finally, Fig. 8 compares the length to width ratio with the experimental data 
for the chick wing from the paper of Haburger and Hamilton [10] (the data 
for the chick leg are similar). Our results are quite consistent with these data. 
We are currently considering experiments which will provide more information 
about the actual time-dependence of the size and shape of the limb [63]. 
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4 Summary 



We have demonstrated that the differential adhesion between cells and dif- 
fusion of morphogens regulating cell growth can generate the paddle shape 
of the chick limb bud. We used a simple model of limb growth with only 
one morphogen (FGF8), and assumed that the cells did not differentiate. The 
next step is to include the condensation of mesenchyme cells into cartilage. 
We know that differential adhesion of cells, together with chemotaxis, leads 
to density- dependent patterns consistent with experiment [31]. Moreover, we 
need to include the anteroposterior asymmetry by introducing the ZPA, and 
use a realistic shape for the AER, with interaction between FGF8 and FGF10 
instead of a constant rate of the FGF8 production. Combining the two mech- 
anisms (condensation and diffusion of morphogens), we hope to obtain a re- 
alistic cartilage pattern in the chick limb (humerus, ulna/radius, and three 
digits). 
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Fig. 3. 3D simulation of the growth of the chick limb bud (2D sections at z = 12). 
The bud consists of: the growing mesenchyme (light blue), AER (yellow), regular 
ectoderm (green), and mesoderm (blue), immersed in medium (white). The width 
of the bud is ~ 4 times smaller than in reality, as expected. 
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Fig. 4. Early stage of 3D limb bud growth (a), and the corresponding FGF8 concen- 
tration distribution (b), shown as 2D sections at z = 12. Only the cells in the apical 
zone with a high concentration of FGF8 grow, generating mainly proximo-distal 
expansion of the bud. 
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Fig. 5. Simulation of excessively fast growth of the chick limb bud. The bud becomes 
bulbous, instead of paddle-like. 
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Fig. 6. Simulation of the growth of the chick limb bud, including the AER growth. 
The shape of the bud is again bulbous. 
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Fig. 7. The time dependence of: the length of the limb bud (a), its width (b), and 
the length of the AER (c). 
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Fig. 8. The time dependence of the length/width ratio (cross signs) compared to 
the experimental values for the chick wing at Hamburger and Hamilton stages 18-23 
(bars). 
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